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ABSTRACT 

In this paper, we develop simple, yet efficient, procedures for sampling approxi- 
mations of the two-Parameter Poisson-Dirichlet Process and the normalized inverse- 
Gaussian process. We compare the efficiency of the new approximations to the corre- 
sponding stick-breaking approximations of the two-parameter Poisson-Dirichlet Pro- 
cess and the normalized inverse-Gaussian process, in which we demonstrate a sub- 
stantial improvement. 
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1 Introduction 

The objective of Bayesian nonparametric inference is to place a prior on the space of prob- 
ability measures. The Dirichlet process, formally introduced in Ferguson (1973), is consid- 
ered the first celebrated example on this space. Ferguson's (1973) original definition of the 
Dirichlet process was based on specifying its finite-dimensional marginals to be Dirichlet 
distributions. An alternative definition of the Dirichlet process, due to Ferguson (1973), 
was relied on normalizing the gamma process. A different constructive definition of the 
Dirichlet process was given by Sethuraman (1994) using a "stick-breaking" approach. We 
refer the reader to the Zarepour and Al Labadi (2012) for more discussion about different 
representations of the Dirichlet process. 

Several alternatives of the Dirichlet process have been proposed in the literature. In this 
paper, we focus on two such priors, namely the two-parameter Poisson-Dirichlet process 
(Pitman and Yor, 1997) and the normalized inverse-Gaussian process (Lijoi, Mena and 
Priinster, 2005). We begin by introducing the stick-breaking definition of the two-parameter 
Poisson-Dirichlet process defined on an arbitrary measurable space (X,A). See for more 
details Pitman and Yor (1997). 

Definition 1. ForO < a < 1, 9 > —a, let (A)i>i be a sequence of independent random 
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variables with a Beta(l — a, 9 + ia) distribution. Define 

i-l 

pi = /9i, j/ i = fc]l(l-p k ),i>2. 

k=l 

Letpi > p 2 > . . . be the ranked values of (p' i ) i>l - Moreover, let (^)»>i be a sequence of 
independent and identically distributed (i.i.d.) random variables with common distribution 
H, independent of (/3j)i>i. Then the random probability measure 

oo 

PbM-) = X>M0 (i-i) 

i=i 

is called a two-parameter Poisson-Dirichlet process on (X, A) with parameters a, 9 and H, 
where 8 X denotes the dirac measure at x. 

It is worth mentioning that Ishwaran and James (200 1 ) referred to the process P' H a 6 (•) = 
Y^H=iP'i$Yi{') as tne Pitman- Yor process, where (p9i>i an( ^ (Yi)i>i are as defined in Def- 
inition [T] The two-parameter Poisson-Dirichlet process with parameters a, 9 and H is 
denoted by PDP(H; a, 9), and we write P H , a ,e ~ PDP(H; a, 9). In the literature, the 
probability measure H is called the base measure of Ph,cx,6> while the parameters a and 9 
are called the discount parameter and the concentration parameter, respectively (Buntine 
and Hutter, 2010; Teh, 2006). The representation ( |1.1[ ) clearly shows that any realization 
of the two-parameter Poisson-Dirichlet process must be a discrete probability measure. 
Note that, the special case PDP(H; 0, 9) represents the Dirichlet process. The law of the 
weights (px,P2, ■ ■ •) is called the two-parameter Poisson-Dirichlet distribution, denoted by 
PD(a, 9). The two-parameter Poisson-Dirichlet distribution has many applications in dif- 
ferent fields such as population genetics, ecology, statistical physics and number theory. 
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See Feng (2010) for more details. On the other hand, the two-parameter Poisson-Dirichlet 
process has been recently used in applications in Bayesian nonparametric statistics such 
as computer science (Teh, 2006), species sampling (Jang, Lee and Lee, 2010; Navarrete, 
Quintana and Miiller, 2008) and genomics (Favaro, Lijoi, Mena and Priinster, 2009). 

The calculations of the moments for the two-parameter Poisson-Dirichlet process are 
carried out in Carleton (1999). Let A be a measurable subset of X. Then 



E(P H , a; e(A)) = H{A) and Var(P H , a>e {A)) = H(A)(1 - H(A))?—^. (1.2) 



It follows from (1.2) that the base measure H plays the role of the center of the process, 
while both a and 9 control the variability of Pn, a ,e around H. Observe that, for any fixed 
set A E A and e > 0, we have 

Pr {\PhMA) ~ H{A)\ >e}< H(A)(1 - H{A))^^- 2 . (1.3) 

Thus, PH, a ,e(A) A H(A) as 9 — > oo (for a fixed a) or as a — > 1, a < 1 (for fixed 6). In 
this paper, "A" and denote convergence in probability and almost sure convergence, 
respectively. 

Analogous to the Dirichlet process, Lijoi, Mena and Priinster (2005) defined the nor- 
malized inverse-Gaussian process P h> q = {P H j(A)} AeA by specifying the distribution 
of (Ph,o(Ai), . . . PH,e(A m )), for a partition A\, . . . , A m of X, as the normalized inverse- 
Gaussian distribution with parameter (9H(Ai), . . . 9H(A m )), where m > 2. See equation 
(3) of Lijoi, Mena and Priinster (2005) for the density of the normalized inverse-Gaussian 
distribution. The normalized inverse-Gaussian process with parameter 9 and H is denoted 
by N-IGP(if; 9), and we write P Hfi ~ N-IGP(if; 9). 
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One of the basic properties of the normalized inverse-Gaussian process is that for any 

AeA, 



E{P Hfi {A)) = H{A) and Var(P H , e (A)) = H{A){1 H{A)) , (1.4) 



where here and throughout this paper 

1 



0VIY-2, 



and T(-2, 6) = J™ r 3 e~*dt 

Observe that, for large 6, £(#) ~ 9 (Abramowitz and Stegun, 1972, Formula 6.5.32, 
page 263), where we use the notation f{9) ~ g(6) if Yuxiq^^ f (9) / g{9) = 1. Like the 



two-parameter Poisson-Dirichlet process, it follows from (1.4) that H plays the role of the 
center of the process, while 9 can be viewed as the concentration parameter. The larger 9 
is, the more likely it is that the realization of Pn,e is close to H. Specifically, for any fixed 

set A G A and e > 0, we have Pn,e(A) A H(A) as 9 — > oo since 

p>1IW-hw>.}< MI d.5) 

Similar to the Dirichlet process, a series representation of the normalized inverse- 
Gaussian process can be easily derived from the Ferguson and Klass representation (1972). 
Specifically, let (EA^ be a sequence of i.i.d. random variables with an exponential distri- 
bution with mean of 1 . Define 



Y i = E 1 + --- + E i . 



(1.6) 
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Let (Yi)i>i be a sequence of i.i.d. random variables with values in X and common dis- 
tribution H, independent of (rV)i>i. Then the normalized inverse-Gaussian process with 
parameter 9 and H can be expressed as a normalized series representation 

PbA-) = E r-uvM ^ (L7) 

j =1 Z^i=l *-* V- i) 



L(x) = I e- t/2 r 3/2 dt, for a; > 0, (1.8) 



X 



where 

L(x) = - 

/2tt 

and 5x denotes the Dirac measure at X (i.e. 5x(B) — 1 if X G B and otherwise). 
Observe that, working with ( |1.7[ ) is difficult in practice because no closed form for the 



inverse of the Levy measure (1.8) exists. Moreover, to determine the random weights in 
( |1.7[ ) an infinite sum must be computed. 

A radically different constructive definition of the normalized inverse-Gaussian process 
was recently established by Favaro, Lijoi and Priinster (2012) using a "stick-breaking" 
approach. Let 1 be i.i.d. random variables with Z± is 1/2-stable random variable with 
scale parameter 1. Define a sequence of dependent random variables (VJ)j>i as follows 

V 1 = V X ' such that X x ~ GIG(0 2 ,l,-i), (1.9) 



X I B 2 1 
VilVx,..., Vi-! = - 1 - such that X t ~ GIG , 1, — ) , i > 2, 

Xi+Zi \ n;=i(i-^) 2 



where the sequences (Xi)i>\ and (^i)t>i are independent and GIG denotes the generalized 
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inverse Gaussian distribution (see equation (2) of Favaro, Lijoi and Priinster, 2012). Define 

s-i 

Pi = V u Pj = V 3 - Vi), j > 2. (1.10) 
i=i 

Moreover, let (li)i>i be a sequence of i.i.d. random variables with common distribution 
H, independent of (Vi)i>i. Define 

oo 

8=1 

Then P H e is a normalized inverse-Gaussian process with parameter 9 and H. 

This paper is organized as follows. In Section 2, we use the stick-breaking represen- 
tations < | 1 . 1| ) and < | 1 - 1 1 p to sample approximations of the two-parameter Poisson-Dirichlet 
Process and the normalized inverse-Gaussian process, respectively. We also show in this 
section that the stick-breaking representations are inefficient for simulation purposes. In 
Sections 3 and 4, we develop simple, yet efficient, algorithms to simulate approximations of 
the two-parameter Poisson-Dirichlet Process and the normalized inverse-Gaussian process, 
respectively. An extensive simulation study evaluating the accuracy of the new methods to 
the stick-breaking approximations is presented in Section 4. The simulation results clearly 
show that the new approximations are more efficient. 



Al Labadi and Zarepour: Interplay of Frequentist and Bayesian 



8 



2 Simulation from Stick-Breaking Representation in Non- 
parametric Bayesian Inference 

Stick-breaking representations are of special interest in Bayesian nonparametric inference. 
For example, they are used in modeling Bayesian hierarchical mixture models (Ishwaran 
and James, 2001; Kottas and Gelfand, 2001). They are also used to compute the moments 
and some theoretical properties of the related priors (Carleton, 1999). In this section, we 
focus on stick-breaking representations of the two-parameter Poisson-Dirichlet process and 
the normalized inverse-Gaussian process from computational point of view. 

The stick-breaking representation of the two-parameter Poisson-Dirichlet process (see 
Definition [T]) can be used to simulate the approximated two-parameter Poisson-Dirichlet 
process using a truncation argument. By truncating the higher order terms in the sum < |1.1[ ), 
we can approximate the stick breaking representation by 

n 
k=l 

where (A)i>i> (jpi)i>x, and (aij)i>i are as given by Definition [T] with f3 n — 1 (hence (5 n does 
not have a beta distribution). The assumption that /3 n — 1 is necessary to make the weights 
add to 1, almost surely (Ishwaran and James, 2001). A random stopping rule for choosing 
n = n(e), where e e (0, 1), is: 

n = inf {i : £ = (1 - ft) ... (1 - p^ft < e} . (2.2) 



The random stoping rule in (2.2) is similar to the one in (2.2) proposed by Muliere and 
Tradella (1998) for the Dirichlet process. The following lemma shows that the weights 
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(pi)i>\ m the stick-breaking representation are not strictly decreasing, almost surely (they 
are only stochastically decreasing). This makes the truncated stick-breaking representation 
inefficient for simulation purposes. 



Lemma 1. Let (p'i)i>i be as in Definition 1 . Then Pr {p' i+ i < p'i} = J Jq f(x, y)dxdy, 
where 

f{x,y) = { ' + X l, x V 1 y \ I{ X > 0}/{0 < y < 1}, 

B(a, b) = r(a)r(6)/r(a + b), a x = 1 - a, p x = 6 + ia and (3 2 = 9 + (1 + 
Proof. Since p\ = A ElfcliC 1 - A), we have 

Pr{p' i+ i<P'i} = Pr{A+i(l-A)<A} 

- p r {,^< 1 }. 

Since pi is a random variable with the Beta(l — a, 9 + ia) distribution, it follows that 
Pi/ (1 — Pi) has the beta distribution of the second kind with parameters a x — 1 — a and 
Pi = 9 + ia (Balakrishnan and Lai, 2009, page 12). That is, /%/ (1 — Pi) has the density 

The lemma follows from the fact that {Pi)i>\ is a sequence of independent random variables 
with a Beta(l — a, # + ia) distribution. □ 

It follows clearly from Lemma [I] that the probability Pr {p' i+ i < p'i] depends on i, a 
and 9. Table 1 depicts some values for this probability. 
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Table 1: Some values of Pr {p' i+1 < Pi}- 



i 


a 


e 


Pr {p' i+l <p'i} 


1 


0.1 


i 


0.672 


10 


0.1 


i 


0.607 


100 


0.1 


i 


0.521 


1 


0.5 


i 


0.598 


10 


0.5 


i 


0.526 


100 


0.5 


i 


0.503 


1 


0.9 


i 


0.5230 


10 


0.9 


i 


0.504 


100 


0.9 


i 


0.500 


1 


0.1 


10 


0.523 


10 


0.1 


10 


0.521 


100 


0.1 


10 


0.511 


1 


0.5 


10 


0.515 


10 


0.5 


10 


0.511 


100 


0.5 


10 


0.503 


1 


0.9 


10 


0.504 


10 


0.9 


10 


0.502 


100 


0.9 


10 


0.500 



Similar to the two-parameter Poisson-Dirichlet process, the stick-breaking representa- 
tion of the normalized inverse- Gaussian process can be used to approximately simulate the 
normalized inverse-Gaussian process process using a truncation argument. By truncating 
the higher order terms in the sum 1[ ), we can approximate the stick breaking represen- 
tation by 

n 

P n ,H,e{-) = I>M-), (2-3) 



k=l 



where (Vi)i>i, (pi)i>i are as given by Definition [I] with V n \Vi,..., V n -i = 1. The assump- 
tion that V n \Vi, . . . , V n -\ = 1 is necessary to make the weights add to 1, almost surely. A 



random stopping rule for choosing n = n(e), where e G (0, 1), is similar to (2.2) with /3, 
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is replaced by V^. Note that, since the (T^)j>i are not independent and the joint density of 
(Vi, Vi + i) is complex for a direct calculation, establishing a lemma similar to Lemma[T]for 
the normalized inverse-Gaussian process is not an easy task. The next table gives values 
of Pr {pi+i < Pi} based on simulation where, the values of the probability is based on 500 
simulated values of each weight with 9 = 1 and n = 50. 

Table 2: Some values of Pr {pi + i < Pi). 



i 


a 


Pr {p i+1 <pi} 


1 


1 


0.536 


10 


1 


0.538 


20 


1 


0.540 


30 


1 


0.548 


40 


1 


0.664 



3 Simulating an Approximation of the Two-Parameter Poisson- 
Dirichlet Process 

The next proposition provides an interesting approach to construct the two-parameter Poisson- 
Dirichlet process. For the proof of the proposition, see Pitman and Yor (1997, Proposition 
22). 

Proposition 1. For < a < 1 and 9 > 0, suppose (pi(0, 6),p 2 {0, 9), . . .) and (pi(a, 0), 
p 2 (a, 0), . . .) has respective distributions PD(0, 8) andPD(a, 0). Independent of (pi(0, 9), 
p 2 (0,0), . . .),let(p\(a,0),p 2 (a,0), . . .) ,i = 1, 2, . . . , be a sequence of independent copies 
of (pi(a, 0),p 2 («, 0), . . .). Let (pi)i>i be the descending order statistics of '{pi(0,9)p t j (a, 0), 
i, j = 1, 2, . . .}. Then (pi,P2, • • •) has a PD(a, 9) distribution. 
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It follows form Proposition [TJ that the weights in the two-parameter Poisson-Dirichlet 
process can be constructed based on two boundary selections of the parameters. The first 
selection is when a = 0. This choice of parameters corresponds to the Dirichlet process. 
The other selection of parameters is when = 0, which yields a measure whose random 
weights are based on a stable law with index < a < 1. Therefore, the Dirichlet process 
Ph,o,o and the stable law process Ph,o,,q are two essential processes in simulating the two- 
parameter Poisson-Dirichlet process. First we consider simulating these two key processes. 

A simple, yet efficient, procedure for approximating the Dirichlet process was was 
recently developed by Zarepour and Al Labadi (2012). Specifically, let X n be a random 
variable with distribution Gamma(6>/n, 1). Define 

r°° i 

G n (x) = Pr(X n >x)= — — e^" 1 ^. t 3 - 1 ) 

Jx T{6/n) 

and 

G-\y) = mf{x:G n (x)>y}. 

Let (Yi)i>i be a sequence of i.i.d. random variables with values in X and common distri- 
bution H, independent of (I"i)j>i. Let 1^ = E\ H V E { , where (i?j)j>i are i.i.d. random 

variables with exponential distribution of mean I, independent of (^)i>i Then as n — > oo, 



i=l 2^=1 G n FT; i=l ^ t=1 V l! 



where N(x) = 9 t l e f dt, x > 0. The right-hand side of (3.2) represents Ferguson's 



representation (1973) of the Dirichlet process. Zarepour and Al Labadi (2012) showed 
that the weights G^ 1 f r^rrj IY^h=i (f^ - ) °^ me new representation given in (3.2) 
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decrease monotonically for any fixed positive integer n. They also provided a strong em- 
pirical evidence that their new representation yields a highly accurate approximation of the 
Dirichlet process. 



The next algorithm uses Zarepour and Al Labadi approximation (3.2) to generate a 
sample from the approximate Dirichlet process with parameters 9 and H . 

Algorithm A: Simulating an approximation of the Dirichlet process. 

(1) Fix a relatively large positive integer n . 

(2) Generate Yi '~ H for % = 1, . . . , n. 

(3) For i — 1, . . . , n + 1, generate Ei from an exponential distribution with mean 1, inde- 
pendent of (V i ) 1 < i < n and let Tj = E x H h E^. 

(4) For i = 1, . . . , n, compute G^ 1 (Ti/T n+ i) , which is simply the quantile function of 
the Gamma(9/n, 1) distribution evaluated at 1 — Ti/T n+1 . 



(5) Set o 6 as defined in (3.2 ) . 



On the other hand, for the stable law process, Pitman and Yor (1997, Proposition 10) 
proved that 

oo —\/ a 

PhM') - E r -i/a M0, (3-3) 

where Ti = E\ + ■ ■ ■ + Ei and (Ei)i>i is a sequence of i.i.d. random variables with an 
exponential distribution with mean of 1. Therefore, the following representation can be 
used to simulate an approximation of the stable law process 

i=l 2^i=l *■ i 
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It is easy to see that the weights [T i 1 ^ a /Yli=i^i ) WG strictly decreasing. 
Thus, simulating the stable law process through the representation ( 3.4[ ) is very efficient. 



The next algorithm can be used to sample from an approximation of the stable law process. 
Algorithm B: Simulating an approximation of the stable law process. 

(1) Fix a relatively large positive integer n. 

(2) Generate Yi '~ H for % = 1, . . . , n. 

(3) For i — 1, . . . , n + 1, generate Ei from an exponential distribution with mean 1, inde- 
pendent of and \etTi = E 1 -\ h 

(4) For each i = 1, . . . , n, the corresponding weights are iy 1 ^ / Ym=i TJ 1 ^- 

Now we present an efficient algorithm for simulating the two-parameter Poisson-Dirichlet 
process. This algorithm is based on Proposition [TJ Algorithm A and Algorithm B. 

Algorithm C: Simulating an approximation of the two-parameter Poisson-Dirichlet 
Process. 

(1) Use Algorithm A to generate n weights of the Dirichlet process. Denote these weights 

by ( Pl (o,e),..., P2 (o,e)). 

(2) Use Algorithm B to generate m weights for an approximation of the stable law process. 
Denote these weights by ( P i(a, 0), . . . ,p m (a, 0)). 

(3) Repeat step (2) to generate n i.i.d. copies of (pi(a, 0), . . . ,p m (a, 0)). Denotes these 
copies by (p\(a, 0), . . . ,j4(a, 0)) , . . . , (p?(a, 0), . . . ,p^(a, 0)). 
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(4) Find the product pi(0, 0)pj(ot, 0) of the weights generated in step (1) and step (3), where 

i = 1, . . . , nand j = 1, . . . , m. That is, find (pi(0, 6)p\{a, 0), . . . ,pi(0, 9)p] n (a 1 0), 

p n (0,e)p n 1 (a,0),...,p n (0,e)pl l (a,0)). 

(5) The weights of the two-parameter Poisson-Dirichlet process are those weights obtained 
in step (4) written in descending order. Denote these weight by (pi)i<i< nm - 

(6) Generate Yj '~' H for i = 1, . . . , nm. 

(7) The approximated two-parameter Poisson-Dirichlet process is given by the representa- 
tion ( |2.1| ) with n in the summation replaced by nm. 

4 Monotonically Decreasing Approximation to the Nor- 
malized Inverse-Gaussian Process 

Mimicking Theorem 1 of Zarepour and Al Labadi (2012), we can construct a similar ap- 
proximation for the normalized inverse-Gaussian process. Specifically, let X n be a random 
variable with distribution IG(0/n, 1). Define 

Q n (x) = Pr(X n > x) = [°° -^r^expt-l (^-+t) + -\ dt. (4.1) 

Jx nV27r L 2 \nH J n J 



Let (li)i>i be a sequence of i.i.d. random variables with values in X and common distri- 
bution H, independent of (r^^i, then as n — > oo 



n Q 1 ( r ' [ oo j, . 
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Here 1^, L(x), and Q n (x), are defined in ( 1.6), ( 1.8 ), and (4.1 ), respectively. 



Remark 1. For any 1 < % < n, rj/r n+1 < r i+1 /r n+1 almost surely. Since Q n 1 is a 
decreasing function, we have Q^ 1 (Vi/V n+ i) > Q^ 1 (Tj + i/T n+ i) almost surely. That is, 



the weights of the new representation (4.2) decrease monotonically for any fixed positive 
integer n. As demonstrated in Zarepour and Al Labadi (2012) for the Dirichlet process, we 
also anticipate that this new representation will yield highly accurate approximations to the 
normalized inverse-Gaussian process. 

Algorithm D: Simulating an approximation of the normalized inverse- Gaussian pro- 
cess. 

(1) Fix a relatively large positive integer n. 

(2) Generate Yj '~' H for i = 1, . . . , n. 

(3) For z = l,...,n + l, generate E { from an exponential distribution with mean 1, inde- 
pendent of (Xi)i<j<„ and let Ti = E 1 H h E*. 

(4) For i = 1, . . . , n, compute Q^ 1 (Ti/T n+ i) , which is simply the quantile function of 
the inverse-Gaussian distribution with parameter a/n and 1 evaluated at 1 — Ti/T n+ i. 
Computing such values is straightforward in R. For example, one may use the package 
' 'GeneralizedHyperbolic" . 
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5 Empirical Results: A Comparison with the Stick-breaking 
Approximation 

In this section, we compare the new approximation of the two-parameter Poisson-Dirichlet 
process (Algorithm C) and the new approximation of the normalized inverse-Gaussian pro- 



cess (Algorithm D) with the corresponding stick-breaking approximations given in (2.1 ) 



and (2.3). First we consider the two-parameter Poisson-Dirichlet process. In the simula- 
tion, we set n = 100, m = 500 in Algorithm C and n = 100 x 500 = 50000 in ( |2~T] ). 
We take H to be the uniform distribution on [0, 1]. We generate 1000 sample paths from 
the two-parameter Poisson-Dirichlet for different values of a and 9 by using the two ap- 
proximations. The sample means of generated processes at x = 0.1, 0.2, . . . , 0.9, 1.0 are 
compared with the true mean of the two-parameter Poisson-Dirichlet process H(x) = x. 
Table 3 shows the maximum mean error (the absolute maximum of the differences between 
the sample means and the true means). For instance, for a = 0.9 and 9 = 10, the maxi- 
mum mean error is 0.00245 in the new approach, while it is 0.01149 in the stick-breaking 
approximation. Similarly, sample standard deviation and the population standard deviation 
can be compared and their maximum errors are reported in Table 3. It is clear from the sim- 
ulation results in Table 3 that both the maximum mean error and the maximum standard 
deviation error in the new approach are smaller than that obtained by the stick-breaking 
approximation. Thus, empirically, simulating the two-parameter Poisson-Dirichlet process 
by using the the new approximation (Algorithm C) gives very accurate results. 
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Table 3: This table reports the maximum mean (max. mean) error and maximum stan- 
dard deviation (max. sd.) error. That is, the absolute maximum difference between 
the approximated mean (standard deviation) and the actual mean (standard deviation) of 

PH,a,e{x) = PH, a ,e((—oo : x\) evaluated at x = 0.1, 0.2, . . . , 0.9, 1.0, where if is a uniform 
distribution on [0, 1]. 









New 




Stick-breaking 


a 


9 


max. mean error max. sd error 


max. mean error max. sd error 


0.1 


1 


0.01155 


0.67082 


0.02332 


0.67082 


0.5 


1 


0.00983 


0.50000 


0.01334 


0.50000 


0.9 


1 


U.UUJ04 


U.ZZoOl 


U.Ul jVj 


U.ZZJOl 


0.1 


10 


0.00993 


0.28604 


0.01097 


0.28604 


0.5 


10 


0.00368 


0.21320 


0.00455 


0.21320 


0.9 


10 


0.00245 


0.09535 


0.01149 


0.13863 


0.1 


50 


0.00236 


0.13284 


0.00341 


0.13284 


0.5 


50 


0.00131 


0.09901 


0.00227 


0.09901 


0.9 


50 


0.00094 


0.04428 


0.00974 


0.05468 
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Figures 1, 2 and 3 show sample paths for the approximate two-parameter Poisson- 
Dirichlet process with a uniform distribution on [0,1] as a base measure with different 
concentration and discount parameters. Distinctly, the new approximation performs very 
well in all cases. On the other hand, a clear disadvantage of the stick-breaking approxima- 
tion appears when a is close to 1 (a < 1). In this case, as seen in Figures 2 and 3, contrary 
to our anticipation (see inequality ( |1.3[ )), the two-parameter Poisson-Dirichlet process is 
not in the proximity of the base measure. Thus, the stick-breaking representation performs 
very poorly when a is close to 1 (a < 1). 

Next we compare the new approximation of the normalized inverse-Gaussian process 
(Algorithm D) with the stick-breaking approximation ( |2.3[ ). In the simulation, we set n = 
50, 9 = 1 and take H to be the uniform distribution on [0, 1]. We generate 500 sample 
paths from the normalized inverse-Gaussian process by using the two approximations. The 
sample means of generated processes at x = 0.1, 0.2, . . . , 0.9, 1.0 are compared with the 
true mean of the normalized inverse-Gaussian process H(x) = v. The maximum mean 
error (the absolute maximum of the differences between the sample means and the true 
means) is 0.01032, while it is 0.01804 in the stick-breaking approximation. Similarly, the 
maximum standard deviation error is 0.06137 in the new approach, while it is 0.07079 
in the stick-breaking approximation. Once again, both the maximum mean error and the 
maximum standard deviation error in the new approach are smaller than those obtained by 
the stick-breaking approximation. 

Figure 4 shows sample paths for the approximate normalized inverse-Gaussian pro- 
cess with a uniform distribution on [0,1] as a base measure with different concentration 
parameters. As seen in this figure, the new approximation performs very well in all case. 



On the other hand, when 9 is large, contrary to our anticipation (see inequality ( 1 .5 1), the 
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Figure 1: Sample paths of the two-parameter Poisson-Dirichlet process Pn,a,e, where H 
is the uniform distribution on [0, 1], 9 = 10 and a = 0.1, 0.5. The solid line denotes the 
cumulative distribution function of H. 



Al Labadi and Zarepour: Interplay of Frequentist and Bayesian 21 



PDP: New 

alpha= 0.8 , theta= 10 



n 
o 



oo 
6 



6 



o 
6 




i — i — i — i — i — i — i — r 

-0.2 0.2 0.6 1.0 



n 
o 



PDP: Stick-breaking 

alpha= 0.8 , theta= 10 



CO 

6 



6 



o 
6 




i — i — i — i — i — i — i — r 

-0.2 0.2 0.6 1.0 



O 



PDP: New 

alpha= 0.99 , theta= 10 



CO 

d 



6 



o 
6 




i i i i i i i r 

-0.2 0.2 0.6 1.0 



n 
o 



PDP: Stick-breaking 

alpha= 0.99 , theta= 10 



CO 

6 



6 



o 
d 




i i i i i i i r 

-0.2 0.2 0.6 1.0 



Figure 2: Sample paths of a two-parameter Poisson-Dirichlet process Pn, a ,e, where H is 
the uniform distribution on [0, 1], 9 = 10 and a = 0.8,0.99. The solid line denotes the 
cumulative distribution function of H. 
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Figure 3: Sample paths of a two-parameter Poisson-Dirichlet process Pn, a ,e, where H is 
the uniform distribution on [0, 1], 9 = 100, 1 and a = 0.99. The solid line denotes the 
cumulative distribution function of H. 
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normalized inverse-Gaussian process is not in the proximity of the base measure. 

Remark 2. When generating samples from the approximate stick-breaking representation 
of the normalized inverse-Gaussian process, the values of the random variables (T^)i>i in 



( 1 .9 ) become numerically 1 after few simulation steps (n < 50). The reason of that, the val- 
ues of the random variables (JQ)i>i becomes very huge comparable to values of the random 
variables {Zi)i>\. In this case, generating samples for the random variable (Xj)j>i from 
the generalized inverse distributions is impossible as these distributions become undefined 
(one of the parameters becomes numerically 0). This problem makes sampling the nor- 
malized inverse-Gaussian process via the stick-breaking representation fails in most cases. 
Our recommendation is to use the stick breaking representation to simulate the normalized 
inverse-Gaussian process only when 6 is very small (9 < 1). 
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